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Abstract. A numerical method is proposed in order to track field lines of three-dimensional divergence 
free fields. Pield lines are computed by a locally valid Hamiltonian mapping, which is computed using 
a symplectic scheme. The method is theoretically valid everywhere but at points where the field is null 
or infinite. Por any three dimensional fiux conservative field for which problematic points are sufficiently 
sparse, a systematic procedure is proposed and implemented. Construction of field lines is achieved by 
means of tracers and the introduction of various Hamiltonians adapted to the "geometrical state" each line 
or tracer is. The states are artificially defined by an a priori given frame of reference and Cartesian coor- 
dinates, and refer to a Hamiltonian which is locally valid at the time step to be computed. This procedure 
ensures the preservation of the volume (fiux condition) during the iteration. This method is first tested 
with an ABC-type flow. Its beneflts when compared to typical Runge-Kutta scheme are demonstrated. 
Potential use of the method to exhibit "coherent" Lagrangian structures in a chaotic setting is shown. An 
illustration to the computation of magnetic fleld lines resulting from a three-dimensional MHD simulation 
is also provided. 

PACS. 5.45.Ac, 05.45.Pq, 47.52.+j, 47.65.Md 



1 Introduction 

With the constant increase of computing power, one is 
now able to simulate more and more precisely complex 
systems, among which three dimensional fluid flows. In 
order to tackle these simulations one is often bound to 
use an Eulerian perspective on flows or flelds. However, 
if one is interested in transport properties, a Lagrangian 
perspective is often best suited For instance the phe- 
nomenon of chaotic advection translates the fact that fluid 
trajectories are chaotic despite the laminar structure of 
the flow j2. In two dimensional flows chaotic advection 
has been studied extensively, as it offers the possibility to 
drastically increase the mixing properties of a given time- 
dependent flow. If the flow is stationnary trajectories of 
passive tracers resume to fleld lines, and no chaos or mix- 
ing occurs in two dimensions. However when considering 
three-dimensional flelds, fleld lines are generically chaotic 
for a stationary flow. Understanding the field Hue chaotic 
structures in these stationary fiows can be therefore con- 
sidered as a first step towards understanding transport 
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in these fiows. Note also that anyhow, a passive tracer 
trajectory is locally tangent to the field fine for the con- 
sidered time, hence computing field lines correctly is as 
well a first numerical step in order to compute trajecto- 
ries. This chaos of field lines has been investigated for 
a long time in plasma physics, especially when concerned 
with the conception of magnetically confining devices such 
as tokamaks 0- Indeed, the conception of such devices 
originates from the fact that charged particles are locally 
trapped by and along magnetic field lines. The chaos of 
field lines is therefore a crucial ingredient for the under- 
standing and characterisation of plasma transport prop- 
erties. Note although that in plasma devices, such that a 
tokamak, numerical studies of three dimensional magnetic 
field lines is simplified thanks to the toroidal geometry and 
a strong anisotropic magnetic field. 

In this paper we propose a volume (fiux) preserving 
numerical approach to the computation of field lines of 
three dimensional divergence free fields, while making no 
assumption on the anisotropy and geometry of the field. 
A first attempt in this direction has been proposed in 0], 
however such schemes are first order ones and are not 
generic to fiux conservative fields. Field lines of conser- 
vative fiux fields as we already mentionned, known to be 
generically chaotic in three dimensions (see for instance 
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(5 for some nice illustrations and references). This prop- 
erty is often illustrated by means of a special local trans- 
formation: the equation of field lines can be written in a 
Hamiltonian form. These Hamiltonians belong to the class 
of 1-^ degrees of freedom systems, which are known to ex- 
hibit generically chaotic behaviour. More specifically the 
field lines exhibit Hamiltonian chaos, implying some vol- 
ume preserving constraints as well as "time-reversibility". 
These features make Hamiltonian chaos quite peculiar, no 
source nor sinks are possibles, neither are attractors, as 
such the dynamics does not have special asymptotics dis- 
playing for instance a fractal attractor. 

Our strategy to compute field lines numerically while 
keeping these specific chaotic features is to use this Hamil- 
tonian mapping and couple it to a symplectic scheme. The 
symplectic scheme ensures the preservation of the fiux of 
the original field, by preserving the volume in phase space 
0. We insist, that this volume preservation is essential 
in Hamiltonian dynamics as it prevents the system from 
having any attractor as well as source or sinks in phase 
space, a feature which is essential when considering trans- 
port properties for large times and computing the kinetics 
limit. In the same spirit, we may expect that paying at- 
tention to volume preservation while computing field lines 
becomes quite relevant when the field itself is "properly" 
computed. We might think therefore to couple this nu- 
meric scheme for field lines to an exactly conservative in- 
tegrator for the Euler fiow, such as the one discussed in 
[3. Finally we would like to mention that the full three 
dimensional incompressible and ideal hydrodynamics can 
be also reduced to a Hamiltonian formaHsm 0, using for 
instance Clebsch variables, but we are then dealing with 
an infinite dimensional phase space, the numerical algo- 
rithms are then not easily implemented although using 
the so called "vortex-line representation" seems promising 
jH] ■ Our approach is much more modest as we do not deal 
with the computation of the field itself, but also allows us 
to consider any divergence free field, such as a magnetic 
field for instance. 

Another potential interest in computing field lines for 
a fiux free field is the possibility to compute Poincare 
sections of the field lines. We are then able to represent 
global information regarding the three dimensional field 
on a plane or a number of planes, offering de facto an- 
other perspective on the fiow and a tool from which phys- 
ical phenomena may have a simpler or different explana- 
tion. In fact, this point of view has been recently used 
quite successfully in the context of identification of three- 
dimensional reconnective structures of a plasma in a pe- 
culiar geometry QHl) which were ill-defined in a Eulerian 
context. With such an approach, we should also be able to 
discriminate the physical importance of islands of regular 
motion within a stochastic sea in the Poincare sections: 
For instance, by identifying the border between chaotic 
and non-chaotic zones, we might identify a localised three 
dimensional structure which does not necessarily have an 
Eulerian counterpart. If the Eulerian counterpart exists, 
we might have a tool to define clearly its shape and we 
might test, as another step, its dynamic robustness. 



The paper is organised as follows. In Section [21 a brief 
summary of field line equations and the Hamiltonian map- 
ping is given. In Sectional we present how the Hamilto- 
nian formaHsm is applied to numerical data. Then in Sec- 
tion 31 we apply the method to two different examples: 
First an ABC like fiow is considered. The field lines are 
computed with a symplectic scheme and the Hamiltonian 
formalism and are compared with a Runge-Kutta inte- 
gration, short-comings of the latter are clearly exhibited. 
A few different choices of parameters are then chosen for 
the ABC fiow, and some three dimensional structures are 
extracted. Finally as a proof of feasibility beyond toy mod- 
els, the method is applied to a magnetic field computed 
directly from an MHD simulation and Poincare sections 
are shown. 



2 Basic Equations 
2.1 Equation of field lines 

Let us briefiy recall the definition of field lines, for this 
purpose let v be a three-dimensional vector field. Field 
lines of v are curves which are tangent to the field at any 
point. This definition may be problematic when the field 
is zero valued at one point, however we may expect for 
a reasonable smooth field that these points if they exist 
may live on a subset of zero measure (their union have 
a zero volume) . In the following we will consider that the 
computed field lines do not cross such a point. In a mathe- 
matical sense we may consider the definitions of field lines 
as: 

V A dM = , (1) 

where dM. stands for a small displacement along the field 
line around a point M. It is usually easier to consider a 
given reference frame and use coordinates, we then rewrite 
Eq.l^l as 

dx dy dz ds , . 

— = = (2) 

Vx Vy V2 V 

where dx,dy, dz are the coordinates of the displacement in 
a given frame, and Vx, Vy, Vz are the coordinates of v(M) 
and V is its norm; ds stands for the norm of dM (we as- 
sumed an orientation of the line is given ) . Let us consider 
now a smooth field, we then choose our coordinate sys- 
tems such that in a given region Vz 7^ 0. Equations ^ are 
directly reduced into a system of two ordinary differential 
equations 

(dx/dz^Vx/vz 

\ dy/dz ^ Vy/vz ' ^ ^ 

where z acts as a time variable. Note that this reduction 
is only locally valid as it may induce some fictitious singu- 
larities if the line ends up in the neighbourhood of points 
where Vz — 0. This point will be discussed again later (see 
sec l3.3ll . Let us now focus on the Hamiltonian formaHsm. 
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2.2 Hamiltonian formalism for divergence free fields 

We are interested in finding a transformation which will 
allow one to describe the evolution of tracers along a field 
line in a Hamiltonian formalism. So let v be a three di- 
mensional divergence free field, such as for instance the 
velocity field of an incompressible fiow or a magnetic field 
and let us reformulate Eq-lEJ into a Hamiltonian formal- 
ism in this frame. The divergence free condition 

V • V = , (4) 
impHes the existence of a potential vector ^ such that 



VAC 



(5) 



C being defined up to a given arbitrary gradient (gauge 
condition) . 

Let us consider the vector potential ^ and relabel its 
coordinates = —p, = H{x, y, z). Now using the gauge 
liberty on the vector potential C, we set — 0. Given 
Eq.® the field v is rewritten as follows 



dH dp dp 
dx dz' dy 



(6) 



Notice taht the equation of motion for the Hamiltonian 
H{x, y, z) are non-canonical. In what follows, we change 
the coordinate system such that they become canonical. 
Given the coordinates {x, y, z), we define q = x and t = z 
and switch to the new system {q, p, r) . We then consider 
the function H{p, q, r) = H{x, y, z). Before moving on, we 
recall that the change of coordinates impHes 



dxf = dqf + dxpdpf 
dyf = v^dpf 



(7) 



where daf = df /da. We then rewrite the equations of the 
field ijSJ in the new coordinate system and readily obtain 
from Eqs.JHl 

, (8) 

[P- dq 

where the dot refers to a derivative with respect to r and 
the tilde was removed from H. The equations are of canon- 
ical Hamiltonian form with the (p, q) acting as canonical 
conjugates, and r acting as a time variable. Hamiltonian 
H is the component of the potential vector associated with 
the chosen fictitious time (the z-direction in our case). 



2.3 Locality 

Note that the singularity for = is not apparent in 
Eq.®. However the singularity remains because the ficti- 
tious time becomes singular for = 0, which if not cared 
about might reverse the arrow of "time". This transfor- 
mation is actually only locally valid. However, first, given 
our hypothesis of a relatively smooth field, it is unlikely 
that we hit a V = point. Second, we have a free choice 



of our coordinate system, thus we can rotate our frame as 
we think is best. With this freedom of transformation, a 
natural choice would be to consider a locally well suited 
coordinate system in order to avoid singularities in Vz - Un- 
fortunately, this strategy is not well suited when comput- 
ing many field lines numerically at the same time. Indeed 
each would need its own local Hamiltonian since the gauge 
condition will be point dependent. And, as a result this 
would be numerically expensive. 

Besides the singularity, injectivity of the transforma- 
tion {x,y,z) {q,p,T) must be ensured if one wants to 
perform time evolution in Hamiltonian variables. Indeed, 
the momentum is formally defined by 



P 



Vzdy , 



(9) 



and, without more specifications, unicity in the determi- 
nation of y given the triplet {q,p,T) is not guaranteed: 
the history of along the field fine is included in p and 
we may recover problems dealt with previously, related to 
Vz = 0. Moreover, it is also obvious in Eq.l|2j that another 
difficulty will arise in a region for which Vy — 0. In this 
region we expect to have dy = and thus a stationary p 
in y, making the transformation likely non-invertible. 

As will be discussed in SecEl we settled for a compro- 
mise between implementability as well as numerical cost 
and debugging time. Let us now discuss one of the key 
reasons for this transformation to be performed. 



2.4 Volume preserving scheme 

We will now use this Hamiltonian structure in order to 
compute numerically the field lines of v, because it will 
allow the use of a symplectic scheme which ensures the 
preservation of invariants, and which given the divergence 
free condition Q , implies the conservation of volume along 
field lines. 

It is well known that intrinsically to any source-free 
systems, there is a volume form of the phase space K", 
say 



/ = dxi A dx2 



A dx„ 



(10) 



such that the dynamic evolution preserves this form. This 
constraints the geometric structure of the field lines. To 
explicit the structure numerically, it is crucial to use a 
scheme preserving this property as shown in section 

For n = 2, source free vector fields are Hamiltonian 
fields and symplectic algorithms are area-preserving maps. 
The situation is different for n > 3: all the conventional 
methods, including the well-known Runge-Kutta and Eu- 
ler methods, are non volume preserving (see |S|llj). By 
means of an "essentially Hamiltonians decomposition" of 
source free vector fields, a volume preserving difference 
scheme was proposed in Ref. [Oj- The algorithm we used - 
we will restrict ourselves to the case n = 3 - does not use 
the decomposition of Ref. : in such a way, computation 
of the "time evolution" of the tracers along the field lines 
requires the explicit introduction of only one Hamiltonian 
(and not two) , the price being the introduction of a time 
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dependent preferential space direction (see Sec. l2.3ll . How- 
ever, the idea is still to use the underlying local Hamilto- 
nian structures, as explained in Sec. l2.2l and integrate the 
dynamical equations by means of symplectic schemes. 
This procedure, indeed, allows the preservation of the vol- 
ume form 

dr 

dx A dy A dz ^ dq A dp A — . (11) 

The meaning of Eq. Ijll|l is the following one : We know 
that dx A dy A dz is preserved as time t evolves. From 
the Hamiltonian framework we have that dq A dp A dr is 
preserved as the fictitious time t evolves which is here 
the coordinate z. The consistency is provided by the way 
we follow the evolutions, i.e., by the nonhnear relation 
dz/dt = Vz between these two evolutions. 



instance, having access to the function H{p, q, r) in or- 
der to compute the Hamiltonian dynamics seems rather 
complex. In order to avoid these difficulties, even though 
the field fines are propagated in the {p, q, r) space, we de- 
cided to compute the —dqH in the {x, y, z) space using 
the expressions 10. 

We therefore have to keep a constant correspondence 
between the two spaces. This is performed by considering 
the expression of p given by Ea. lll4|l . and explains our 
gauge choice to keep it as simple as possible. In order to get 
(x, y, z) given a value of (p, q, r), only y is given implicitly 
by Ea. ljl4ll since x — q and z — t. The inversion is done 
numericafiy using a simple Newton scheme, for which the 
first guess for the root is given by the Euler method. 



3.3 Dealing with pseudo-singularities 



3 Application to numerical data 
3.1 Fourier representation 

We now present a possible way to implement this for- 
malism when dealing with a field obtained by numerical 
simulation. We consider a quite common case of a field 
with periodic boundary conditions. If this is the case, it 
is a common practice to consider the Fourier space for 
the computation as for instance when one uses a pseudo- 
spectral method. We therefore consider a case where all 
the Fourier components Vk of v are known. Here the vector 
k is used as a mode index. We shall also for convenience 
consider that the number of modes is finite, which is more 
realistic when dealing with numerical data. Using equa- 
tions (EJ we then write the Hamiltonian and momentum 
in the {x, y, z) coordinates: 



H{x,y, 



fc„#o 



ik-x 



E 

1 A;^=0 



E 



P = 



fey#0 



2^ ^Ky^ 



„ik 



iky 



, fe„=0 



Vk.zB 



Note that the condition = is not fully sufficient to de- 
termine H and p, and the equations and I114II which 
are defined up to a function of x and z, refiect our choice 
of trying to keep p as simple as possible, leaving the "com- 
plexity" in H. 



3.2 Numerical scheme 

In order to compute field lines, we take advantage of the 
Hamiltonian formalism and use a symplectic scheme. In 
the present case we consider the fourth and sixth order 
Gauss-Legendre scheme proposed in ^2] • However the trans- 
formation from the {x, y, z) to (p, g, r) is not trivial. For 



Up to now we have assumed that the Hamiltonian trans- 
formation is globally valid. By the choice of coordinate 
system, the direction z acts as a time variable, and is 
therefore assumed to be a monotonic function along the 
field line. One may therefore wonder what happens when 
we hit a fold, which turns the field line in the reversed 
direction. If this is the case z can not be monotonic and 
the Hamiltonian transformation fails. Our way of dealing 
with such problems is to assign a state to each field fine, 
corresponding to which direction the fine is propagated 
in. Indeed when we reach a point where Vz = 0, one may 
expect that Vx or Vy is not null. Since our coordinate sys- 
tem is arbitrary, we may conceptually rotate our frame in 
order to get Vz ^ m this new coordinate system and 
propagate our field with the symplectic scheme. Ideally, 
it would be best suited to locally rotate the coordinate 
system and try to take the arc-length along the fine as 
a measure of the fictitious time variable. However when 
(X2)computing many field lines together, this procedure would 
mean to keep track of all coordinate systems and rotate 
them after each step to the local Frenet system before per- 
\ forming the transformation and the next time step. This 
Cl'|)seems numerically expensive. We settled for a more prag- 
/ matic point of view. Field lines are thus assigned a number 
i € {1, 2, 3} corresponding to which coordinate system the 
^■j^^^Hamiltonian transformation is appHed to. These numbers 
corresponding to a cyclic permutation of the coordinates 
{x, y, z) of a coordinate system given once and for all. 
Furthermore as was previously discussed, problems may 



also arise when Vy = 0. If such a seldom situation occurs, 
we simply permute the current x and y directions. For- 
mally this corresponds to a different choice for p and q and 
a different choice of gauge for ^. A new number j G {1, 2} 
is therefore assigned to the field line. 

When changing the coordinate system, we have to make 
sure that we will continue to "advance" along the line, 
meaning we have to be careful and make sure that the 
new Hamiltonian transformation will not end up in going 
backtrack along the already computed path. We therefore 
choose an orientation for the fines and impose the condi- 
tion 

V ■ dM > . (15) 
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It is obviously important to ensure that this condition is 
still satisfied when changing coordinates. Since it is con- 
ceptually easier to deal with a positive time step dr, and 
since field lines are invariant when multiplied by a scalar 
field, we perform the change v — v when necessary in 
order to comply with Eg . 111511 . This implies that another 
number k G {1, 2} must be assigned to the field line. 

At last, since we try to avoid performing a rotation 
of our coordinate system, except for simple permutation 
of the coordinates, we may hit points for which the field 
is along one of the coordinates, leaving us with two zero- 
valued coordinates and problems of invertibility. In order 
to deal with such situations, we actually have to consider 
a second reference frame, not issued by a permutation. 
We may then consider for instance a it /A rotation around 
the z-direction of the original coordinate system, and deal 
with problematic points in this frame. Given all possibili- 
ties while computing the field line, its state (the numbers 
(i, j, k)) can thus take twelve different values, and we can 
be in two different coordinate systems, giving us a broad 
range of twenty four Hamiltonians in order to compute a 
field line. 



4 Applications 

In order to test the algorithm and start to explore the 
possibilities offered by the analysis of field line chaos, we 
started with a basic three dimensional fiow, namely the 
ABC-fiow p3.14 . We first validate our numerical scheme 
by considering a choice of parameters for which no singu- 
larity occurs and the Hamiltonian transformation is global. 
Further on, we show the advantages of using a symplectic 
scheme over a non-symplectic one. We then consider inte- 
grable and chaotic fiows for which singularities occur. In 
the latter case, Lagrangian structures are exhibited as a re- 
sult of the analysis of Poincare sections. Finally as a proof 
of concept we compute the magnetic field fines, in the con- 
text of the dynamo process, of a magneto-hydrodynamic 
fiow. 



4.1 ABC flow 

In order to test the algorithm we consider the following 
periodical fiow which may be seen as a particular example 
of the well know ABC fiow [TMl] : 

Vx = COS y — e sin z 

Vy = sinx + ecos.z (16) 
Vz = COS X — sin y -t- wqz , 

where e and wqz are parameters. 

We started to benchmark the code by comparing its re- 
sults with the one obtained using a non symplectic fourth 
order Runge-Kutta scheme on Equations (pj. Note that 
in order to have some consistency, comparison is done 
with the symplectic fourth-order Gauss-Legendre integra- 
tor 22]. For this first test we took e = 0.9 and uoz = 4 in 



order to avoid singularities which occurs at points where 
Vz — 0. Since Wz > 0, z is used as the time variable, and 
due to the periodic character of the fiow, we can reduce 
our analysis of the fiow on the Poincare map defined by the 
crossing of field lines with the successive z = [27r] planes. 
The resulting plot is depicted in Fig.Qlusing the two meth- 
ods. The "time step" used for the run was 5z = tt/IOO. To 
the naked eye, the same section is obtained for both cases, 
points solely differ in the chaotic region. The chaotic be- 
haviour of steady field fines in three dimensions is clear. 
This first test is used in order to confirm that the proposed 
method computes what it is expected to. 

Now let us consider the advantages of using such a 
symplectic algorithm. For this purpose we consider the 
integrable case e = 0, and we compare the evolution of 
a field line starting from the same initial condition (for 
visual purposes) but computed using both methods. Since 
the fiow is integrable we expect to see a closed fine on the 
Poincare map. 

First we consider both cases with an identical rela- 
tively large time step 5z = tt/IO. The results are shown 
in Fig. 12 and clearly demonstrate that the non symplec- 
tic method has some shortcomings as the expected closed 
line becomes a cloud of points, and inward diffusion is 
observed. Now in order to be "fair", since the non sym- 
plectic scheme is almost twice faster on an identical com- 
puter than the "conserving" one, we took a smaller time 
step 5z = 7r/20 for the Runge-Kutta method in order 
to match the integration times (CPU times) between the 
methods. Still, as seen in the insert in Fig. 12 some in- 
ward diffusion is observed with the Runge-Kutta scheme. 
These non-conservative features may induce spurious ef- 
fects, for instance when considering transport properties 
in the kinetic fimit, while transport may be not Gaussian 
(see for instance jlSiTZj) one may end observing a diffu- 
sive behaviour with a diffusion coefficient governed by the 
algorithm, as this diffusive behaviour wifi end up eras- 
ing in a finite time any potential memory effect generated 
by the presence of regular trajectories in a mixed phase 
space. In this sense using a conservative algorithm despite 
its extra numerical overhead, we may avoid false interpre- 
tations which may arise with Runge-Kutta as for instance 
possible transient anomalous behaviour. 

In the rest of the paper we use the sixth-order Gauss- 
Legendre scheme. 

4.2 Flow with pseudo singularities 

We now test the code on a fiow for which the problems of 
singularities and non invertibility occur. Namely we con- 
sider a series of fiows for which vqz 0. We first con- 
sider the integrable case (e = 0). The results are shown 
in Fig. El It is easy to see from Eas. (jl6|l that the fiow be- 
comes somewhat two-dimensional, as the z dependence of 
the fiow completely vanishes. For this "two-dimensional" 
steady fiow cos x+siny = H2D is a constant along the field 
line and corresponds to the stream function which acts as 
a Hamiltonian. This property gives rise to the closed lines 
observed in the Poincare section depicted in Fig. El Thus 
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for any value vqz > we expect to have identical sec- 
tions, which is the case as we take different values of vqz 
ranging from 3 to 0.01 {vqz 0) observed in Fig. El An 
additional feature occurs when vqz — 0, all field lines ap- 
pear to become periodic and localised (forming a closed 
line). In this regard, we like to pinpoint a specific fact 
about our choice for the equations of the ABC flow given 
by Eas. ljl6ll . In fact it is usually a common procedure to 
take Vz — H2D(x,y) (see for instance |15ll8j V In the in- 
tegrable situation this procedure leads to a constant Vz, 
and helical field lines. In these regards, the choice of the 
minus sign in the third equation of Ens. ((T6ll drastically 
changes the behaviour of field lines. From the perspective 
of testing the numerical scheme, this choice also proved to 
be more challenging than usually JH], since the field lines 
"bite their tails". 

Since the numerics are giving us what we expect, we 
may now move to a non-integrable situation. For this pur- 
pose we consider e — 0.15 with a vqz — 0: the integrable 
case with closed field lines is perturbed. A Poincare section 
is shown in Fig. 31 One can observe regions with regular 
quasi-periodic field lines, and regions of chaos. In the reg- 
ular region, one can notice some delocalised field lines, 
as well as apparently localised ones (at least in the z- 
direction). In the chaotic sea, one notices a region with a 
vanishing density of points, in other words a region which 
is not often crossed. Since the depicted section is com- 
puted in the x = [27r] plane, one can expect that regions 
in these planes for which Vx — Q are not crossed. It implies 
that regions around the curves given by cos y = £ sin z are 
not often crossed. These curves correspond to the light re- 
gion in Fig. ^ for one, and for the "centre" of the regular 
region in Fig. 31 for the other one. Before concluding on 
this particular example, we want to mention the existence 
of two islands in the chaotic region nearby the low-density 
region that the method is able to detect. 



4.3 Lagrangian structures 

In fact the Poincare section shown in Fig.01is quite generic 
for systems with one and a half degrees of freedom, with 
regions of regular motion and a chaotic sea. In this section 
we describe the shape of the field lines corresponding to 
the regular regions. Indeed contrarily to a true phase space 
section from a given Hamiltonian, the section depicted in 
Fig. 31 corresponds to field lines obtained typically by a 
"product" of local Hamiltonians. The "time" dependence 
is therefore not always monotonic. We thus may expect 
more complexity than the typical braid of helical struc- 
tures which is expected in the phase-space extended with 
a time direction. Besides this possible increase of complex- 
ity, the regular structures may have also a crucial physical 
impact. For instance, if one consider the field v as being 
the velocity field of an incompressible fiuid, the transport 
properties of passive scalars will be along the field lines. 
Thus the properties of these regular regions may greatly 
infiuence global properties of transport or mixing in this 
fiow. 



In order to localise the coherent Lagrangian structures 
associated with regular field lines, we just consider initial 
conditions within a regular zone in the Poincare section, 
and computed the resulting field lines. To our knowledge 
this type of methodology has not been applied to generic 
divergence free fields, in particular to the considered ABC 
fiow. Results are shown in Fig. Typically one notices 
two types of behaviours, localised field lines, which are 
drawing a finite surface, and unlocalised ones. Note the 
mushroom-like shape of the unlocalised field line, which 
if considered as the trajectory of a passive tracer, may 
appear as successive sheets in an experiment. Among the 
unlocalised regular field lines, one can also find a sim- 
ple helical like form, when one considers an initial con- 
dition within the two small islands located not far from 
the no crossing region. To conclude on these structures, 
it is worthwhile mentioning that even by looking at the 
sections trough different planes {y = Q [27r] or z = [27r]) 
the same structures are obtained. They just have differ- 
ent sections with these planes. One may also infer the 
Poincare plot depicted in Fig. that the localised struc- 
tures (Fig. 13 lower panel) appear as sandwiched between 
the unlocalised mushroom ones (Fig. upper panel). We 
may also note that the observed "coherent structures" do 
not have necessarily a simple tubular shape, hence we may 
speculate that using this technique to uncover spatially ex- 
tended coherent structures in more complex fiows such as 
MHD, may yields some unexpected shapes with maybe im- 
portant physical consequences and selection mechanisms. 

4.4 Magnetic Field from MHD 

We shall now consider a full feature numerical study and 
move on to the magnetic field obtained through magneto- 
hydrodynamics (MHD) simulations. The origin of cosmic 
magnetic fields can be investigated in the frame of MHD. 
The problem is to extract the underlying mechanisms of 
the dynamo process, i.e. of the self-excitation of a mag- 
netic field by plasma or fiuid motions. To describe the 
generation of large-scale magnetic fields from small-scale 
turbulence many mechanisms are involved and, in the fol- 
lowing, we will focus on the so-called a^-dynamo (see for 
instance |20]1. 

The equations of incompressible MHD, written in the 
usual units, read 

D , 
— u = -Vp + b • Vb + i/V^u + f 

-^b = b • Vu + ?/V^b (17) 

Vb = Vu = 0, 

where by definition, -q is the magnetic diffusivity, v the 
kinematic viscosity and the scalar p is the sum of the hy- 
drostatic and the magnetic pressures. We choose v = -q = 
. Here, both fiuid velocity u and magnetic field b are 
expressed in Alfven speed units. The fiow is forced by a 
divergence free field at scales much smaller than the width 
L = 27r of the cubic and periodic box. a^-dynamo is en- 
hanced by imposing the forcing to be maximally helical: 
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its relative helicity is one at any time jl9l2f)| . Initial con- 
ditions are u = and b ~ 0. The flow is driven by the 
forcing and it turns out that the Reynolds number based 
on the Taylor microscale is around a few units: this is a 
moderate turbulent regime. In fact, a noteworthy effect of 
the Lorentz force is to decrease by roughly a factor ten 
the Reynolds value. In Fig. are drawn time evolutions 
of the total kinetic and magnetic energies of the plasma. 
Clearly, owing to the kinetic forcing, kinetic energy is ini- 
tially rapidly accumulated. The growth of the magnetic 
energy initiates then a decrease of the kinetic energy and 
a quasi-saturated state is reached at time T = 300. The 
eddy turnover time is typically of a few time units. Trans- 
fer of helicity from the forcing to the flow and from the 
flow to the magnetic fleld, but also Alfven effects and an 
inverse cascade process (after time T ~ 15 here) lead to 
the generation of a robust anisotropic large scale magnetic 
fleld, characterised by a pile up of magnetic energy at the 
largest scales Pl 20J. Once the large scale magnetic fleld 
reaches a stationary state, a rough snapshot of the Eule- 
rian fleld can be summarised by the approximation 

b cos zGx -I- sin zey + e bk exp(ik • x) , (18) 

k|~/c_p 

where e is a small parameter to indicate that energetically 
dominant modes are at the largest scales k — 2t: / L = 1 
and fc_F is the forcing scale. In Fig.d we draw a Poincare 
section of the three-dimensional magnetic fleld lines result- 
ing from the magnetic fleld obtained at time T — 900. This 
is the section of 100 magnetic fleld lines. Initial positions of 
the tracers to determine the flelds lines are randomly dis- 
tributed in the periodic box. The computation is done by 
including in the fleld an energetically consequent fraction 
of the modes. Inclusion of all the modes would have been 
too computer time consuming and, since in this paper we 
put our effort on showing a proof of concept and we do 
not focus on the details of the Lagrangian structures. For 
instance, dissipative Eulerian structures are not included 
in the structure of the fleld. The existence of chaotic La- 
grangian structures is however clearly observed in Figd 
It is also observed a"drift of the structures" according to 
their height z. In fact, this drift trivially originates from 
the Eulerian fleld as can be seen by inspection of the flrst 
term on the r.h.s of Eq. Ijl8|l . We like to insist on the fact 
that Fig. 13 clearly suggests the existence of a non trivial 
underlying Lagrangian and chaotic organisation in duality 
with the well established coherent Eulerian structures sug- 
gested by the approximation ifTsji . These remarks suggest 
that the study of the Lagrangian structures in this type 
of flow offers promising prospects. Indeed, chaos of fleld 
lines of a flow, neither necessarily turbulent nor unstation- 
nary, is an ingredient known to clearly enhance a dynamo 
effect; an clear example being the G.O. Roberts dynamo 
where an ABC flow is used[22]. The large scale term of 
the saturated magnetic fleld in the simulations is, in fact, 
also of ABC type. However the role of the turbulence in 
the generation of the former, but also how the inverse cas- 
cade generates helical structures in a turbulent flow are 
not clear. Thinking those processes in terms of Lagrangian 



structures in parallel with a classiflcation and measure of 
their helicities might allow a new understanding of the 
interplay of turbulent, helical and chaotic phenomena in 
dynamo processes. Finally, we should emphasise that the 
presence of chaos does not necessarily rely onto any com- 
plexity in the small scale perturbative term of Eq. Ijl8|l . 
Indeed, it probably essentially results in the presence of 
the large scale Beltrami structure, bg=o = 'tVAbg^o with 
K = — 1 =Cte, as a consequence of the Arnold theorem|13j. 
We may also emphasise on the fact that the existence of 
chaotic fleld lines is in turn a signature of the three dimen- 
sional nature of the fleld, and considering the anti-dynamo 
theorems in two-dimensional flows we may expect that the 
mechanisms underlying the dynamo effects are generically 
correlated with the existence of chaotic fleld lines 



5 Conclusion 



The main purpose of this work is to propose a way to com- 
pute fleld lines of divergence free flelds, while keeping the 
conservative flux condition valid. Field lines are computed 
by performing a locally valid Hamiltonian transformation, 
which in turn is numerically computed using a symplec- 
tic scheme. The method is theoretically valid anywhere 
but at points where the fleld is null or inflnite. For flelds 
for which such points are sufficiently sparse that we con- 
sider they numerically do not exist, a systematic proce- 
dure has been proposed. Field lines are propagated using 
up to twenty four different Hamiltonians, depending in 
which "state" each line is. The states being artiflcially de- 
flned by an a priori given frame of reference and Cartesian 
coordinates, and refers to which Hamiltonian is currently 
locally valid for the time step to be computed. The poten- 
tial application to three dimensional flelds with periodic 
boundary conditions is discussed and explicit usable ex- 
pressions for the Hamiltonian and momentum are given. 
We have also provided a detailed possible actual potential 
implementation into a numerical code of the ideas, and 
presented numerical results obtained for speciflc systems. 
In this procedure, we have shown some advantages of us- 
ing this Hamiltonian-symplectic formalism over a regular 
Runge-Kutta scheme for the fleld lines. 

Moreover, while testing the numerical scheme, fleld 
lines of ABC type flows have been investigated. For a con- 
sidered integrable flow with no mean flow, we have shown 
that fleld lines can be all closed. Also a Poincare section 
of a non-integrable flow, for which the fleld lines can not 
be described by a global Hamiltonian has been computed 
and a mixed phase space with a chaotic sea and regular re- 
gions is exhibited. From this analysis Lagrangian coherent 
structures are extracted, and either localised or extended 
structures have been found. Finally in order to show an 
application of the numerical tool in order to tackle physi- 
cal problems, the fleld lines of the magnetic fleld resulting 
from an magneto-hydrodynamic simulations in a dynamo 
process have been computed. 
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Fig. 1. Poincare section of the flow l|16|l with voz ~ 4 and 
e = 0.9. The Poincare section is performed on the z — [27r] 
plane. The flow is computed using the 4th-order symplectic 
scheme with a time step 5z — tt/IOO. 
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Fig. 2. Poincare section of the flow l|16|l with «oz = 3 and 
e = 0. The Poincare section is performed on the z — Q [27r] 
plane for only one initial condition (3.142,4.742) ~ (7r,37r/2). 
The black line corresponds to the expected quasi-periodic tra- 
jectory obtained with the fourth-order symplectic scheme. The 
cloud of green dots corresponds to the trajectory obtained by 
the fourth-order Runge-Kutta scheme, showing some diffusion 
towards the centre of the island. Integration is performed up 
to z = TT 10^, with a time step Sz = vr/lO. In the insert the 
same simulation is performed for the Runge-Kutta scheme but 
with a smaller time step 5z = 7r/20 diffusion is still present, 
the dark line being a replot of the results obtained with the 
symplectic scheme and the larger time step 5z — tt/IO. 
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Fig. 3. Left: Poincare section of the flow llfill for the integrable 
case e = 0: note that the same section is obtained as we take 
different values of uoz ranging from 3 to 0.01 {voz — > 0). The 
* signs, correspond to points obtained for vqz = 0, the fleld 
lines appear to be periodic for this case. The Poincare section 
is performed on the z — [2tt] plane with St — 7r/50. Right: 
Pour fleld lines computed for the integrable case £ = and 
voz = 0. Pield lines are indeed periodic and form closed lines. 
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Fig. 5. Upper panel: regular drifting fleld line. Lower panel: 
Localised fleld line. 



Fig. 4. Poincare section of the flow Hl(i|l for a chaotic case 
e — 0.15, for a flow with pseudo-singularities. One can infer 
a chaotic region and a regular one. One also notices a region 
which is not often crossed in the chaotic sea (lighter zone). The 
Poincare section is performed on the a; = [27r] plane. The plot 
is made with 100 fleld lines which are iterated 2 10^ times with 
a time step of &z = vr/lOO. 
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Fig. 6. Time evolution of kinetic energy (green line) and mag- 
netic energy (blue line). 
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